/* 
 * File:   main.cpp
 * Author: Bangonkali
 *
 * Created on March 23, 2012, 8:14 PM
 */
#include <iostream>
#include <cmath>
#include <cstdlib>

// based on discussions from here http://j.mp/GILnX3
// example taken from here http://j.mp/GILphB

using namespace std;

int gauss_jordan(int N, double A[3][4], double r[3]) {
    double scale;
    
    for (int p = 0; p < N; p++) {
        double divider = A[p][p];
        if ( divider != 0 )
        {
            for (int i = 0; i <= N; i++) {
                if (i != p) {
                    scale = A[i][p];
                    for (int j = 0; j <= N; j++) {
                        A[i][j] = (A[p][j] * scale / divider) - A[i][j];
    //                    printMatrix (3, A);
                    }
                }
            }
        }
    }
    
    for (int p = 0; p < N; p++) {
        scale = A[p][p];
        for (int j = 0; j <= N; j++) {
//            cout << "A[p][j] = ";
//            cout << "A[" << p << "][" << j << "] / ";
//            cout << "A[" << p << "][" << p << "] = ";
//            cout << A[p][j] << "/" << scale << " = ";

            if (A[p][j] != 0 && scale != 0)
                A[p][j] = A[p][j] / scale;
            else
                A[p][j] = 0;

//            cout << A[p][j];
//            cout << endl;
//            printMatrix (3, A);
            
        }
        
        r[p] = A[p][N];
    }
    return 0;
}

//int main(int argc, char** argv) {
//    double A[20][21]; // input
//    double X[20]; // result
//    
//    A[0][0] = 1; 
//    A[0][1] = 2; 
//    A[0][2] = 3; 
//    A[0][3] = 3;
//
//    A[1][0] = -3; 
//    A[1][1] = 1; 
//    A[1][2] = 5; 
//    A[1][3] = -2; 
//
//    A[2][0] = 2; 
//    A[2][1] = 4; 
//    A[2][2] = -1; 
//    A[2][3] = -1; 
//    
//    printMatrix (3, A);
//    gauss_jordan(3, A, X);
//    printMatrix (3, A);
//    
//    for (int i=0; i<3; i++) cout << X[i] << "  ";
//    return 0;
//}
